Criticality in transport through the quantum Ising chain 
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We consider thermal transport between two reservoirs coupled by a quantum Ising chain as a 
model for non-equilibrium physics induced in quantum-critical many-body systems. By deriving 
rate equations based on exact expressions for the quasiparticle pairs generated during the transport, 
we observe signatures of the underlying quantum phase transition in the steady-state energy current 
already at finite reservoir temperatures. 
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Quantum phase transitions (QPTs) are drastic mani- 
festations of quantum fluctuations in many-body systems 
that lead to critical behavior and states of matter with 
symmetry-broken phases [1] . Traditionally, such changes 
of state are associated with changes in ground state prop- 
erties at zero temperature. Recently, considerable atten- 
tion has turned to the role of excited states, as these 
become important if the fate of QPTs under modifica- 
tions such as finite temperature or non-equilibrium due 
to dissipation, external control, or driving is addressed. 
Among the fundamental issues then is the very definition 
of a QPT under non-equilibrium conditions, and further- 
more the question of how non-analyticities of the ground 
state, phase diagrams, critical exponents, scaling behav- 
ior, or the order of the transition are modified. 

From the point of view of recent advances in simula- 
tions and realizations of QPTs in, e.g., cold atoms [2, 3], 
these are urgent issues, as experimental systems always 
contain a certain amount of dissipation and require ex- 
ternal control and read-out mechanisms. For exam- 
ple, a classical external control allows one to explore 
novel states of matter and effective interactions which 
are absent in equilibrium [4-7] and lead to new types 
of phase diagrams with additional phases, transitions 
lines and critical points. Similar results have been ob- 
tained for QPTs in open, dissipative systems, either de- 
scribed by additional dissipative terms in equations of 
motions (similar to mean- field- type laser equations [8, 9]) 
or Lindblad-type master equations, e.g. [10]. A further 
line of theoretical work concerning excited-state phase 
transitions has also mostly been restricted to mean-field 
type QPTs (such as the Lipkin-Meshkov-Glick [11] or the 
Dicke-superradiance [12] model), where non-analyticities 
in quantities like the density of states or excited state en- 
ergies are essentially related to singularities in a (classi- 
cal) energy landscape. At zero system temperature, these 
can be made visible by suitable measurement devices - 
as has e.g. been suggested for the Dicke model [13]. 

In this letter, we investigate whether such features are 
also accessible in the non-equilibrium regime at the ex- 
ample of the quantum Ising chain in a transverse mag- 
netic field. Its interesting properties can either be stud- 
ied in existing Ising ferromagnets [14] or in corresponding 



quantum simulators, e.g. based on NMR techniques [15], 
trapped ions [16], atomic spins [17, 18], or electrons float- 
ing on liquid helium [19]. Quite some work has been de- 
voted to thermal transport along open spin chains [21- 
24]. Complementing these approaches, we consider here a 
scenario where transport occurs perpendicularly through 
a closed ring of spins. Going beyond linear response [25], 
this enables us to explore the extreme non-equilibrium 
regime. We use the chain, which has a transition from a 
paramagnetic into a Zi symmetry-broken ferromagnetic 
phase with long-range magnetic ordering at zero tem- 
perature, as a thermal coupling between two reservoirs 
whose temperature gradient drives a stationary, thermal 
current. Combining exact results for the Ising model 
with the weak-coupling, Born-Markov-Secular type de- 
scription of open systems [20] , we demonstrate that the 
thermal current displays criticality at the same value of 
the Ising chain coupling parameter as for the zero tem- 
perature (ground state) QPT. This has to be contrasted 
with stationary equilibrium quantities such as the aver- 
age energy density or magnetization that display critical 
behavior only at T = 0. 

Traditionally, the vicinity of the critical point of equi- 
librium QPTs at finite temperatures has been character- 
ized by spatio-temporal order parameter correlation func- 
tions, which in particular for the quantum Ising chain are 
known to display a relatively complex behavior [1]. We 
find that the stationary non-equilibrium current through 
the critical system itself may serve as a powerful tool in 
order to reveal quantum critical properties (which docs of 
course not exclude the possibility of much more complex 
features to be found in, e.g., temporal non-equilibrium 
correlations such as current noise spectra). 

We also find that neither the position nor the order 
of the phase transition arc changed by adding dissipa- 
tion, provided the Lindblad operators are evaluated at 
the correct frequencies, i.e., in positivity-conserving secu- 
lar approximation. Our conclusion from this observation 
is that phenomenological extensions of QPTs which sim- 
ply add constant dissipation terms on top of equilibrium 
models have to be treated with caution. 

Model. — Our full model Hamiltonian rl = T~Ls +"Hsb + 
Hb consists of the quantum Ising chain in a transverse 
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field for N spins 



N 



AT 



(1) 



where 3 describes the coupling to an external magnetic 
field, J the inter-chain coupling to nearest neighbors, and 
periodic boundary conditions are assumed <J Z N+1 = erf. 
The model is analytically diagonalizable for finite N and 
displays a second order quantum phase transition be- 
tween a paramagnetic phase (for g > J) and a ferromag- 
netic one (g < J) [26]. Furthermore, transport through 
the chain is enabled by coupling to two equilibrium reser- 
voirs Hb = J2a=s,DJ2 q UJ aqbl iq b aq , source (S) and drain 
(D), at different temperatures. Here, we consider non- 
interacting bosons (e.g., photons or phonons) with mo- 
mentum q and energy ui aq (we set H = 1), having in mind 
a thermal-transport type flow of energy that is mediated 
by boson-induced spin-flips via the interaction 



H S B = a ? i h U b U + i&fi"*] 



(2) 



with interaction constants h % aq that depend on the spe- 
cific realization of the model in an experimental con- 
text. We mention that the derivations below can also 
be carried out for fermionic reservoirs albeit with some 
modifications due to finite chemical potentials that oc- 
cur in that case. More importantly, we choose a par- 
ticularly simple case of "Hsb that is amenable for ana- 
lytical calculation by assuming the coupling strengths as 
site- independent h l aq — > h aq . Effectively, the coupling 
to the Ising chain is then via the large spin operator 

We first diagonalize 'Hs in the usual way, apply- 
ing a Jordan- Wigner, Fourier, and Bogoliubov trans- 
formation ([26, 27]) which, in the subspace of an even 
number of quasi-particles and for even N (these con- 
straints will be tacitly assumed further-on), transforms 
the system Hamiltonian into = J^ fe tk{l\lk — 1/2) 
with fermionic annihilation operators "f k [27]. The 
quasi-momentum k may assume half-integer values 
k = ±1/2, ±3/2, . . . , ±(N - l)/2. The single-particle en- 
ergies are defined by 



e k = 2fi^ (1 - S f + s 2 - 2s(l - s) cos , (3) 

where we have introduced a dimensionlcss phase param- 
eter by fixing Qs = J and Q(l — s) = g with energy scale 
0, see Fig. 1. 

The very same transformation maps the interaction 
Hamiltonian, Eq. (2), to [27] 

^sb = j/Vl-2^ M 2 1 + (K| 2 - M 2 ) lllk (4) 

^ k - 

+u k v„ k (ll k ^_ k + 7-fc7+fe) > [ h «i b U + h - c -] ' 

- J rv a 
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FIG. 1. (Color Online) Spectrum of the quantum Ising chain 
of length TV = 6 versus phase parameter s = J/Q — 1 — g/Q 
emphasizing the subspace of even quasi-particle numbers and 
vanishing pair momenta. In the thermodynamic limit TV — > 
00, the second order quantum phase transition comes with 
a vanishing energy gap above the ground state, cf. Eq. (3). 
Possible transitions between states for second order perturba- 
tion in h ak are marked by vertical lines with corresponding 
single-pair energies 2et . Thin dotted curves in the background 
complete the full spectrum of Eq. (1). 



where the coefficients are defined by v k oc s sin (^S^) and 
u k oc (l — s — scos (^jir-) + efc/(2f2)) with normalization 

\u k \ 2 + \Vk\ = 1. 

Obviously, the (spatially homogeneous) interaction 
"Hsb does not couple subspaces with different values of 
the total momentum P — J^ fe kjt'yk, since it either does 
not create particles at all (first line in Eq. (4)) or creates 
or annihilates only quasiparticle pairs of opposite mo- 
menta. In the paramagnetic phase with s <C s C rit, these 
are spin waves running in opposite directions that can 
be interpreted as two hard-core bosons created or anni- 
hilated by the absorption or emission of a single reservoir 
photon/phonon. In the ferromagnetic domain s 3> s cr it, 
the quasi-particles correspond to kinks between domains 
of opposite spin directions. 

Assuming that the system is initially prepared in 
the subspace of vanishing total momentum (e.g., in its 
ground state |0)), it now suffices to consider the subspace 
of pairs with opposite quasi-momenta only. In this sub- 
space, the basis elements can be conveniently constructed 
from the ground state via 
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Til , T13 , 



n 

fc>0 



10) ,(5) 



where n k e {0, 1} denotes the occupation of a 
quasi-particle pair with momenta (+k,—k) such that 

(7fc7fe + 7l fc 7-fc) I") = 2n fc \n). 

Rate Equation. — Applying lowest order perturbation 
theory in the coupling strength and in the relevant sub- 
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space (employing Born, Markov, and secular approxima- 
tions [20] in the standard way appropriate for open sys- 
tems) yields a rate equation 



Pn = X] I S ^nm Pm 

m \ a / 



(0) 



for populations of the system energy cigenstates 
Pn = {n\ P where the transition rates Inm 

due to reservoir a admit only creation or anni- 
hilation of single quasi-particle pairs, see vertical 
lines in Fig. 1. Assuming thermal reservoir states, 
the transition rates (n ^ m) evaluate to [27] 
Inm = r«(Amn) [1 + n a (A mn )] \{n\ M x \m)\ 2 with 
energy differences Amn = Em — En, system energies 
En = J2k>o e fc(2«.fc — 1), spectral coupling densities 
r a (uj) = 27r^ 9 \h aq \ 2 5(uj — uj q ), and bosonic occupations 

n a (uj) = ^e l3aUJ — l] 1 (diagonal values j nn follow from 
trace conservation). 

When both reservoirs are held at sufficiently low tem- 
peratures such that n a (2e k ) ~ for all k > 3/2, 
the system will relax to the subspace spanned by the 
ground state |0) = |0, 0, . . . , 0) and the first excited state 
|1) = |1, 0, . . . , 0). All states of higher occupations relax 
via the annihilation of higher-momentum quasi-particle 
pairs towards this subspace, see Fig. 1. The reduced tran- 
sition rates are proportional to the matrix element [27] 



i<oiM,|i>r 



,s 2 sin 2 



N i 



1-2 S (1- S ) [1 + cob(#)] 



(7) 



with its second derivative with respect to s becoming 
nonanalytic at s cr ;t = 1/2 as N — > oo. In the continuum 
limit however, the dynamics cannot be constrained to 
the 2x2 subspace anymore as many levels approach the 
ground state, cf. Eq. (3), which motivates study of the 
full-dimensional problem. 

Using Eq. (6) and the rates j nm , we obtain an ana- 
lytical result for the steady state solution [27], 



Pn = H 

fc>0 



[n(2e fc )]" fc [l + n(2 £fc )] ] 
l + 2n(2e k ) 



(8) 



which - similar to Rcf. [28] - is completely governed 
by an effective average bosonic occupation n(u>) = 



— j2 — r — . However, our system has more than one 
allowed transition frequency, which implies that the sta- 
tionary state (8) is non-thermal (i.e., cannot be described 
by a single effective temperature) as soon as the reservoir 
temperatures arc different (rig (a;) =^ nu(uj)). Eq. (8) en- 
ables us to calculate the stationary values of the energy, 
the magnetization, and the current both for finite chain 
lengths and in the continuum limit N — > 00. 

Energy and Magnetization. — For the mean stationary 



energy, wc find [27] 



1/2 



E 



fe>0 



N-yoo 



2n(2e k ) 



N 



l + 2n(2e(«)) 



(Lk,(9) 



where we have introduced the continuum of system ener- 
gies c(k) = E(jvk)' At zero temperature, where n(2e(«;)) = 
0, the system settles to the ground state, and the en- 
ergy density can be expressed by a complete elliptic inte- 
gral E/N -> — 2p£e(4s(l - s)), with a divergent second 
derivative at s cr it = 1/2. This divergence, which reflects 
the usual ground state QPT criticality of the Ising chain, 
is also predictable from analyzing the analytic structure 
of the integrand in (9) at zero temperature. For finite 
temperature and also in non-equilibrium setups where 
n(u) > for all u> > 0, the energy density and its higher 
derivatives remain analytic. 

We find a similar behavior for the magnetization, 
which for large N becomes [27] (v(k) = i'(jvfc)) 



N 



1/2 



1 - 4 



\v(k)\ 2 + n{2e{n)) 
1 + 2n(2e(/c)) 



dn 



(10) 



At zero temperature, the integral is similarly solved by 
normal elliptic integrals and those of the first kind, which 
display a divergence in the first derivative of the magne- 
tization density with respect to s. However, at finite 
temperature all derivatives of the magnetization density 
are analytic, which is most evident in the trivial high- 
temperature case where n{2e(n)) — > 00. 

Heat Current. — This changes drastically, however, 
when we consider the heat current through the Ising 
chain from one reservoir to the other. Analysis of the 
transition rates (e.g., by introducing energy counting 
fields as in [29]) yields our main result for the current 
of net emitted bosons at the drain, 



I = Y.( Em - E n)lnmPm (11) 
n, m 

= \ - 2e k A 2 k T s (2e k )T D (2e k ) [n s (2e k ) - n D (2e k )} 
f^ Q T s {2e k ) [1 + 2n s (2e fc )] + r D (2e fe ) [1 + 2n D {2e k )\ ' 

where the second line follows after a straightforward cal- 
culation by inserting the stationary state and explicitly 
evaluating the transition rates [27]. Here, we have intro- 
duced A k = Au k v k = sin (^). 

Evidently, the current is antisymmetric when S O D, 
vanishes at equilibrium, and is positive when the source 
temperature exceeds the drain temperature (which im- 
plies ns(oj) > ni)(w)). Most important however, in the 
thermodynamic limit TV — > 00 the current I directly re- 
flects the signatures of the ground state quantum phase 
transition of the Ising chain. Formally, this correspon- 
dence is visible by the integral representation of /, which 
shows a divergence of its second derivative with respect 
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FIG. 2. (Color Online) Renormalized energy current / and 
its second derivative w.r.t. s (inset) versus phase parameter s 
for different chain lengths N = 4, 40, oo (dotted, dashed, and 
bold solid, respectively) and for different source temperatures 
fi/3s = 0.1, 0.5, 1.0 (black/brown, red/orange, and dark/light 
green, respectively). The dash-dotted purple curve denotes 
the analytically accessible case of noiy) — > 0, ns(u) — ¥ oo, 
and N — >• oo. Other parameters: Q/3d = 10, Ts(efc) = 

r D (e k ) = r. 



to the phase parameter s at all temperatures, see Fig. 2. 
The second derivative of the integrand in the continuum 

1/2 

version I/N = L j(s, k)cIk of Eq. (11) will at the crit- 



ical point s cr ;t 
d 2 j{s,K) 



d 2 i 



1/2 for small k diverge as [27] 

32nr s T D (p D -p s ) 



,1/2 



tt{TsI3d + T d Ps)k 



whilst the integrand itself and its first derivative remain 
finite. Even for the extreme non-equilibrium, infinite 
thermobias regime the analytically obtainable [27] cur- 
rent displays a divergence of the second derivative at 
Scrit = 1/2, compare the dash-dotted curves in Fig. 2. 

Conclusion. — For transport through a closed Ising 
chain homogeneously coupled to two bosonic thermal 
reservoirs we have analytically shown that signatures 
of the underlying QPT are manifest in the energy cur- 
rent already at finite temperature and also in the non- 
equilibrium regime - in contrast to system observables 
as mean energy or mean magnetization. For finite sizes 
N, all quantities remain analytic but precursors of the 
QPT are already visible at rather moderate chain lengths. 
Slightly perturbing the coupling symmetry, all subspaces 
of the model - cf. thin curves in Fig. 1 - may be weakly 
coupled, but since for near-multistablc systems the sep- 
arate current contributions are additive [30], we expect 
the sensitivity of the current to the underlying QPT to 
remain. Furthermore, in the weak-coupling limit con- 
sidered here, the position of the phase transition is nei- 
ther changed nor are new phases introduced by the cou- 
pling to reservoirs. We note that higher order transport 



cumulants (noise or time-dependent current autocorrela- 
tion functions) can be expected to yield even more com- 
plex features and thus might prove useful tools to analyze 
QPTs out of equilibrium. 
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Supplementary Material 



i 



Equations in the main manuscript are referenced by square brackets. 

Mapping the system Hamiltonian to non-interacting fermions 

The Jordan- Wigner transform (JWT) 



n x — 1 - 2<-t c 



<=-( Cn+c t) n (1-2 



(i) 



m— 1 



maps the spin-1/2 Pauli matrices non-locally to fermionic operators c m . Inserting the JWT into the Hamiltonian [1] 
in the main manuscript, one has to treat the boundary term with special care 



JV-l 



n s = -.g53( 1 ~ 2c « c ") ~ J 13 ( c ™ + c «)( c »+i + c l+i)( 1 - 2 4c«) - J(c N + c^ N ) 



n=l 
N 



n=l 
JV-l 



'N — l 



I (1 - 2cJ, c„) 



(ci+4) 



-9^](l - 2c+c„) - J 5^(4 - c «)(4+i + c n +i) + J(c]y - C/v)(4 + Ci) 



iV 



1(1 -2ctc„) 



(2) 



where we have extensively used the fermionic anticommutation relations. Introducing the projection operators on the 
subspaces with even (+) and odd (-) total number of fcrmion quasiparticlcs 



1 



N 



i ± n a 2 C t c „) 



(3) 



we may also write the Hamiltonian (2) W$ = + V )'Hs{'P + +V ). It is straightforward to see that terms with 
different projectors and with n < N vanish 



= V+(l - 2cic n )V- = P"(l - 24 c n )V+ , 

= 7>+(4 - c„)(4 +1 + c n+1 )V- = P"(4 - c„)(4 +1 + c n+1 )V+ 

For the boundary terms one finds similarly 

N 



(4) 



(P++p-)(4-c iV )(4 + c 1 ) 



1(1-240 



(73- 



= + p-)(4 - c Jv)(4 + Cl )(2T + - 1)(V+ + V~) 
= P+(4 - cw)(c t 1 + Cl )V+ - p-(4f - c N )(c\ + Cl )V- , 



(5) 



such that we can finally write the Hamiltonian (2) as the sum of two non-interacting parts with either an even or an 
odd total number of fermionic quasiparticles 



V 



N 



TV 



g^i 1 ~ 2 4c„) - J 53(4 - C„)(4+l + C n+l) 
n—l n—l 

N N 

-9 53( 1 ~ 2c « c ») ~ J 53 ( c « ~ c «)(4+i + c n+i) 



v- 



n=l 



n=l 



(6) 



Note that this requires to define antiperiodic boundary conditions in the even (+) subspace 4^ = — and periodic 

boundary conditions in the odd (-) subspace 4+i = + c i ■ 

Since the even subspace is relevant to our model, we further seek to diagonalize the Hamiltonian 



N 



N 



-.9 53( 1 ~ 2c n C n) ~ J 53 ( C « ~ C «)( C «+1 + C «+l- 



(7) 



n=l 



2 



with antiperiodic boundary conditions cjv+i = —c\. Translational invariancc suggests to use the discrete Fourier 
transform (DFT, preserving the anticommutation relations due to its unitarity by construction) 



„ — I7T/4 



N 



(8) 



which is compatible with the antiperiodic boundary conditions when N is even such that k takes half-integer values 

fce{±i ±l±\,...}, where \k\ < ^ (9) 
The DFT maps the Hamiltonian into 



Hg = -gNl + {% ~ Jcos(k2Tr/N)]5ld k + J sin(k2n/N) 



(10) 



Now, the observation that only positive and negative frequencies couple (conservation of one-dimensional quasi- 
momcntum), suggests to use the reduced Bogoliubov transform 



c k = u +kl+k + v*_ k l-k , 



(11) 



which mixes positive and negative momenta and where the a priori unknown coefficients have already been labeled 
suggestively (a more general ansatz would eventually of course yield the same solution). Since the new operators ^ k 
should be fermionic, we obtain from the orthonormality conditions 



1 



\u+ k \ 2 + \v- k \ 2 , = u +k v* +k + U- k v*_ k = {v* +k ,v*_ k ) ( j 



(12) 



Demanding that the Bogoliubov transform eliminates all non-diagonal terms (of the form "f- k j+ k etc.) yields (by 
combining positive and negative k) the equation 



g- J cos fc- 



CI = 2 

= (v_ fe ,u_ fc ) 



N 



(u +k v- k - U-kV+k) + 2J sin ( fe^ ) (u- k u +k + v- k v +k ) 



-2 [g- Jcos(fe^)] +2Jsin(fe^) 



-2Jsin(fc^) 



-2[g- Jcos(fcf )] J \v+ 



u +k 



(v- k ,u-k)M 



u +k 

V+k 



(13) 



All equations can be fulfilled when we choose (u+ k , v +k ) T as the normalized positive energy eigenstate of the matrix 
M. with eigenvalue 



and {v*_ k ,u*_ k Y 



e+ = +2y / g 2 + J 2 - 2gJcos(k2ir/N) = e k 
(— v* +kl +u* +k ) T as its negative energy eigenstate with eigenvalue e k 



(14) 



-2a/ g 2 + J 2 — 2g J cos(k2ir / N) . To be more explicit, we have 



u k = 



v k = 



g - Jcos(k2n/N) + yjg 2 + J 2 - 2gJ cos(k2n/N) 

g - Jcos{k2n/N) + V ' g 2 + J 2 - 2g J cos{k2i: / N) + [Jsin(fc27r/iV)] 2 

Jsin(fc27r/iV) 



(15) 



g - Jcos{k2n/N) + g 2 + J 2 - 2gJcos(k2ir/N) + [Jsin(fc27r/iV)] 2 
Using these solutions, we obtain when N is even 



*J = E W + J2 - 2 sJcos (fcf) ( 7 t 7fe - = E 



(16) 



which reproduces the inline equation before [3] in the main manuscript. Note that the case of odd N and/or diago- 
nalization in the odd subspace may be treated similarly. 
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Mapping of the interaction Hamiltonian 

Obviously, the used transformations do not affect the reservoir part, such that it suffices to transform M x = X)i=i a f 
with the very same transformations as before. Inserting the JWT (1) yields 

JV 

M x = Nl - 2 J2 cUn ■ (17) 

n=l 

Furthermore, inserting the DFT (8) leads to 

M x =Nl-2^2c{c k . (18) 
fe 

Finally, inserting the Bogoliubov transformation (11), replacing k — > — k in some terms and exploiting that the 
coefficients (15) are real yields 



M x =Nl-2 

k 



E [\ u k\ 2 "fhk + \v k \ 2 lkll +u k v^ k (j+kl-k + 7-fe7+fc) , (19) 



which by using the fcrmionic anticommutation relations is equivalent to Eq. [4] in the main manuscript. For later 
convenience we write this as a sum over positive k- values only 

A4=iVl-4^| Ufc | 2 l-2^(l-2|« fc | 2 ) (lh k +"fi k "f- k ) +4^ MfcVfc (- 7 | fc7 l fe +7 +fe7 _ fe ) . (20) 

k >0 k>0 k>0 



Derivation of the rate equation 

We rely on previous results in the literature that yield for an interaction Hamiltonian of the form 'Hsb = A<8> B 
under Born, Markov, and secular approximations a Lindblad-type master equation. When the spectrum of the 
system Hamiltonian is non-degenerate (and here more specifically, when states coupled in the master equation are 
energetically non-degenerate) , this Lindblad master equation couples only the diagonal elements of the density matrix 
in the system energy eigenbasis to each other, i.e., it can be written in the form of a rate equation 



Pa 



^ Jab,abPbb — ^ Jba,baPaa , (21) 



where a, b label the different system energy eigenstatcs. Note that the refined condition of non-degeneracy is for finite 
N always fulfilled, as e.g. the intersection point in Fig. 1 between 1 001) and 1 1 10) in the main manuscript are between 
uncoupled states. The transition rates are given by 

lab,ab = l(E b -E a )\(a\A\b)\ 2 , (22) 

where 7 (w) is the Fourier transform of the bath correlation function, 7 (w) = J C(T)e +1LJT dr 

I (e + ' MsT Be~' MsT B) e +iUT dr. Specifically for our model there is only one, which for two thermal source and drain 

reservoirs becomes 

C (T)= E ([h aq bi q e + ^+h.c] [h a , q 'bl q ,+h.c])=J2J2\ h <*«\ 2 [( b U^ 

aa'qq' a q 

oo 

= ^ E / r <» KM e+i " r + (1 + n«M)e-^ r ] du; , (23) 
n a o 

where we have introduced for ui > the spectral coupling density T a (ui) = 27r^ 9 \h aq \ 2 S(uj — io aq ) (recall that 

u a q > 0) and n a (uj) = \e^ aUJ — l] 1 denotes the Bose distribution for reservoir a held at inverse temperature /3 a . 
Analytically continuing the spectral coupling densities to negative frequencies via T a (— u) = — r a (+o;) and exploiting 
that n a {— lu) = —(1 + n a (+uj)) yields after a simple integral transformation 

c ( T ) = ^E/ r »[ 1+n »] e+i " Tdtu . ( 24 ) 



4 



which enables to directly read off the Fourier transform 7(0;) = J2 a F a {u) [1 + n a (u)] = S a 7a(k>). Evidently, the 
contributions of the two reservoirs enter additively in the rate equations, such that by labeling the energy eigenstatcs 
in the relevant subspace (a —> n) we recover the rate equation with its coefficients stated below Eq. [6] in the main 
manuscript. 



Low Temperature Limit 



At sufficiently low temperatures, such that n a (2ek) <C 1 for all k > 3/2 but still n a (2e 1 / 2 ) = 0(1), it is evident 
from Fig. 1 that most excited states will relax towards the two lowest states |0) = |0, . . . 0) and |1, . . . 0) = |1). The 
dynamics in this subspace is governed by the rate equation 



Po 
Pi 

with the matrix element 



= |(0|M ;c |l)| 2 ^r Q (2 ei/2 ) 



|<0|M x |l)| a =4u? /2 u? /2 



-n Q (2e 1/2 ) 
+n Q (2e 1/2 ) 



1 + n a (2e 1/2 ) 
1 + n a (2e 1/2 ) 



Po 
Pi 



S 2 sin 2 (£) 



1/2 1/2 ~ 1 - 2,(1 - a) [l + cos(f)]' 
compare Eq. [7] in the main text. Consequently, the current in this effective low-temperature limit becomes 

I 



T s (2ey 2 )T D (2e 1/2 ) [n s (2e 1/2 ) - n D (2e 1/2 )] 



s 2 sin 2 



(f) 



r S (2ci/ 2 ) [1 + 2n s (2e 1/2 )] + r D (2e 1/2 ) [l + 2n D (2e 1/2 )] 1 - 2s(l - a) [l + cos (f )] 



(25) 



(26) 



(27) 



which modifies the usual bosonic current through a two-level system by the matrix element |(0| M x |1)| . Eventually, 
the s-dcpcndcncc of this prefactor yields the non-monotonous dependence of the current on the phase-parameter s in 
the low-tcmpcraturc curves in Fig. 2. 



Non-equilibrium Stationary State 

The stationary solution of the rate equation can even for non-equilibrium (different temperature) configurations be 
obtained using basically two ingredients. First, we note that the Fourier transforms of the bath correlation functions 
obey the usual Kubo-Martin-Schwinger conditions j a (—uj) = e~^ aU} j a (+oj), which lead when the system is coupled 
to only one bath (e.g. by setting Td(lo) = 0) to thcrmalization of the system with the temperature of the remaining 
reservoir (e.g. Pg 1 )- Formally, such a thermal state is characterized by the ratio of diagonal elements to be 

£H. = e -«%-%! = »(£n - Em) ^ , 2g , 
Pm 1 + n(E n - Em) ' 

where n(w) corresponds to the Bosc distribution of the connected reservoir. For coupling to multiple reservoirs we 
use that the occupations of the different reservoirs enter linearly and just weighted by the different tunneling rates to 
motivate the ansatz (A nm = En — Em) 

pn n(A n m) . T s (uj)ns(uj) + T D (u))n D {uj) 

— = : : -, A r , n[u) = . (29) 

Pm l + n(A nm ) r s (w) +T d (lj) 

Indeed, one can easily prove for the rate equation 

Pn= J2 J2 r ^ Amn ^ 1 + n ^ Amn ^^ n ^ M ^ m ^ 2 P m 

5Zr Q (A nm )[l + n ct (A nm )]|(m|A4|n)| 2 | p n (•'!()) 



1 m^n a 

the validity of the stationary state by inserting 



Pm — : rr-r zPn — ^ „ rrr- tt rrPn yoi-) 

l + n(A m n) L a T a{A m n) [1 + n a (A m n)\ 
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and using that r o (A m n) = -r a (Anm) and n a (A mn ) = — [1 + n a (A nrn )]. By a sequence of pair annihilations - 
compare Fig. 1 in the main manuscript - it therefore follows that any stationary occupation may be connected to the 
ground state occupation po via 

- - tt ( n(2e k ) 
Pn = Po[[[ T 



fe>0 



+ n(2e fe ) 



(32) 



The latter is fixed by the normalization Tr {pn} = 1 



e - e n(i 

"1/2=0 n(jv_i)/2=0 fc>0 

which yields Eq. [8] in the manuscript. 



n{2e k ) 
+ n(2e k ) 



k>0 Ln k =0 



n(2e k ) 
+ n(2e k 



Po 



fe>0 



1 + 2n{2e k ) 
n(2e k ) 



(33) 



Stationary Energy 



Rewriting the system Hamiltonian (16) as 



E e * hhk + 7i fe 7-fe - 1 



fc>0 



(34) 



implies for its diagonal matrix elements in the relevant subspace (n| "Hg \n) = ^2 k> o e k(2n k — 1). The stationary 
expectation value of the system energy then becomes with Eq. [8] in the main text 



n 



(E) = Tr {H+p} = Y, <™l ?4 I") Pn = E ^ E( 2n * ~ *W 

k>o n 

I 1— 71 k 

■ (2n fc -l) = E r 



^ ^ [n{2e k )] nk [l + n(2e fc ) 



1 + 2n(2e fc ) 



fc>0 



2n(2e ft ) 



(35) 



where we have used that X^ fe =o " fc 7+2^. — ~ = ^ holds for each k separately in the second line. In the thermodynamic 
limit (TV — > oo) and noting that all relevant quantities actually depend onit = k/N, the sum is easily converted into 
an integral, and we recover Eq. [9] in the main text. 



Stationary Magnetization 



Similarly, we evaluate the diagonal matrix elements of the magnetization operator (20) 

(n\ M x \n) = TV - 4^ \v k \ 2 - 4^T (K| 2 - | Ufe | 2 ) n k = TV - 4^ [|« fe | 2 + (l - 2|w fe | 2 ) n k 



k>0 k>0 

which can be inserted in the stationary expectation value to yield 



fe>0 



(M x ) = i n \ Af» I") Pn = N - 4 



n 



fe>0 



^~4^|^| 2 -4^(l-2 



k>0 



fe>0 



4^(l-2| Wfc | 2 ) ]T n >< 

k>0 n k =0 

n{2e k ) 



[n(2e k )] nk [l + n(2e fc )] ] 
1 + 2n(2e fe ) 



1 + 2n(2e fc ) 



AT A V N + n ( 2g fc) 

A? l + 2n(2e fc ) 



(36) 



(37) 



Finally, the sum over k can similarly be converted into an integral to yield Eq. [10] in the main text. Further- 
more, by inserting the coefficient (15) in the continuum representation and zero-temperature limit, we obtain for the 
magnetization density 



1/2 

(m x ) = =1-4 J v 2 ( K )dn 

o 



gg(4g(l - a)) + (1 - 2 5 )£ g (4s(l - a)) 

7r(l - S) 



(38) 



where £e(x) and Sk{x) denote the complete elliptic integral and the complete elliptic integral of the first kind, 
respectively. 
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Stationary Current 

The stationary current of bosons emitted to the drain can for example be obtained by inserting energy counting 
fields in the off-diagonal matrix elements of the rate equation matrix, i.e., to perform in Eq. (30) the replacements 

rz3(A mn ) [1 + n D (Amn)} -> r D (A mn ) [1 + n D (A m n)} e +iA mnx ? (39) 

which automatically takes into account that Amn > corresponds to emission into the drain and Amn < to 
absorption. Note that in the latter case one would use Tjj(— x) [1 + Ud(— #)] = T d{+x)tid{+x) ■ This upgrades the 
rate equation by a counting field p — C(x)Pi an d the stationary current can then be obtained from the stationary 
state by deriving the rate matrix with respect to the counting field x 

1= (-i)Tr {£'(())/>} =J2 ^mnr D {Amn)[l + n D (Amn)}\(n\M x \m)\ 2 pm 
n m^n 

Arnn^D{Arnn)[l + n D (Amn)}\(n\M x \m)\ 2 pm 

nm ■. A mn >o 

Y A nm T D (A n m)n D (Anm)\(n\ M x \m)\ 2 p m 
nm : A nm >o 

= \2e k m k T D (2e k ) [1 + n D (2e k )] (Au k v k ) 2 p m ~ 2e fe (l - m k )T D (2c k )n D (2e k ) (Au k v k f p m 

m fc>o 



^ 2e k T D (2e k )(Au k v k ) 2 2J [m k [1 + n D (2e k )} - (1 - m k )n D (2e k )] pm 



k>o m 

1 



£ 2e k T D (2e k )(4u k v k ) 2 ]T [m k [1 + n D (2e k )} - (1 - m k )n D (2e k )\ M 2e ^ " t 



1 + 2h(2e k ) 

k>0 m k =0 v K/ 



J2^kr D (2e k )(4u k v k ) : 



n(2e k ) - n D (2e k ) 



1 + 2n(2e fe ) 



_V- 9 u s 2 r s (2e k )T D (2e k ) [n s (2e k ) - n D (2e k )} 

~ j^ Q Zek[mkVk) T s (2e k ) [1 + 2n s (2e k )} + T D (2e k ) [1 + 2n D {2e k )] ' [W) 

which with evaluating the prcfactor A k = Au k v k from (15) becomes Eq. [11] in the main manuscript. 
The continuum representation of the current becomes (in wide-band limit T a = T a (2e k )) 

1 _oo ^ Wsin 2 ^™) T s T D [ns(2e( K ))~n D (2e(K))] _ T w 

N ~ J e(n) T s [l + 2n s (2e( K ))]+T D [l + 2n D (2e(K))] dh - J ^ K > dK - ^> 



At the critical point and for small k, the integrand behaves as 



9 7 ^ 



--1/2 



32*tl(p D -0 s )r D T s 2 

„ „ k + C{k}, (42) 

1 sPd + 1 dps 



which together with Eq. [12] in the main manuscript leads to divergence of the second derivative of the current at the 
critical point for all temperature configurations. 

This can also be seen in closed form in the infinite thermobias regime (ns{2e(n)) — > 00 and no{2t{n)) — > 0), 
where (41) becomes 



1/2 

I , ^0 f sin 2 (27TK) , 



T [(1 - 2s(l - a))£ E (4a(l - a)) - (1 - 2s) 2 £ K (4 s (l - a))] , (43) 
where £e{x) represents the complete elliptic integral and £k(%) the complete elliptic integral of the first kind. 



